NPS ARCHIVE 
1969 

DUBAC, C. 

DERIVATION THE OPTIMAL CONTROL FOR AN 
ALL-WEATHER AIRPLANE LANDING SYSTEM 

by 



Carl Hugo Dubac 



‘A *N - 

jUGNIfl J13HS 

mr°f*5 



DUDLEY KNOX LIBRARY 
NAVAL POSTGRADUATE SCHOOL 
MONTEREY, CA 93943-5101 



United States 

Naval Postgraduate School 




THESIS 



DERIVATION OF THE OPTIMAL CONTROL 
FOR AN ALL-WEATHER AIRPLANE LANDING SYSTEM 



by 



Carl Hugo Dubac 



June 1969 



Tkii document kcu> been approved f$oA pubtic \e- 
ieaAe and baJLe; -it6 di6Vu.bation -u untimLted. 




library 

NAVAL POSTG RADUATE school 

M-JN t’lklY. calif. 93940 



Derivation of the Optimal Control 
For An All-Weather Airplane Landing System 



by 



Carl Hugo JDubac 

Major, United States^Marine Corps 
B.S., University of Michigan, 1956 



Submitted in partial fulfillment of the 
requirements for the degree of 



MASTER OF SCIENCE IN ENGINEERING ELECTRONICS 

from the 

NAVAL POSTGRADUATE SCHOOL 
June 1969 



z 



i r ^ n \ ' 




•uiaqsAs aqq jo aoueurtogjed 
aqq e^GJ^snxiT qoxqM pa^uasaxd axe saxduiBxg ’AxoqoaCGjq 
aqeqs paxxsap e Buxuxjap Aq pue 'uoxqxpuoo qqbxxj uinxxq 
-Txxnba ue ^noqe paziieauTt st qoxqM auBxdJXB aqq xog x s P OUI 
XeoxqBUiaqqBui e fiuxdoxaAap Aq uiaxqoxd BuTqaexq leauix g sb 
pa^Bxnuijog: st buxpuBX eqi, *auBXct:rTB rfr-d SBXbnoa pxeuoaow 
aqq ut fiuxpuBX Jaqqea/A~xx e ue J° sascqd x eu TJ sq3 J°J 
TaXX o:c: luoo b aATjap oq pasn st Axoaqq x o:i: ) uod Xeuixqdo 



iDVHiSav 



DUDLEY KNOX LIBRARY 
NAVAL POSTGRADUATE SCHOOL 
MONTEREY, CA 93943*5101 



Thesis by: Carl Hugo Dubac entitled Derivation of the Optimal 

Control for an All-Weather Airplane Landing System. 

ERRATA 



Page Line 
21 



Change 

On Fig. 4 change q to Q = q 



24 



8 



Insert after Eq. (3.3): where 







x a 


= x , ,v , 

03 O 


x * = 


X <^o' Z « = 


Z, > V 
CO 1 






z a 


= Z • v , 

03 O 


= 


M, v , and 
03 o 


M* : 
a 


26 


8 


Change X , 


x « 


' Z q' 


to 








x , 

q 


x a 


' X 6 e ' V 


• • • 


35 


13 


Eq. 


(4.11) 


should read 





x(t) = Ax(t) + B_u(t) + c . 



rr^qasjxon* Y3 n 1 
' :HC? :T V > ' 

I t . .V 1 













TCPA'J JATE SCHOOL 

w./l; ,'l* v . calif 



TABLE OF CONTENTS 



I . INTRODUCTION 9 

II. PROBLEM FORMULATION 11 

III. PLANT MODELING 19 

IV. OPTIMAL CONTROL THEORY 30 

V. SPECIFICATIONS, DESIRED TRAJECTORIES, 38 

PERFORMANCE MEASURE 

VI. INVESTIGATION PROCEDURES 57 

VII. RESULTS 63 

VIII. CONCLUSIONS 89 

APPENDIX A. REPRESENTATIVE DATA FOR THE F-4J 91 



AIRPLANE IN THE LANDING 
CONFIGURATION 

APPENDIX B. SUMMARY OF NUMERICAL VALUES 92 
COMPUTER PROGRAM 93 
LIST OF REFERENCES 102 
INITIAL DISTRIBUTION LIST 103 
FORM DD 1473 105 



3 



LIST OF ILLUSTRATIONS 



Figure Page 

1. Category II ILS Beam Accuracy Requirements 13 

2. Phases of Automatic Landing 16 

3. Airplane in Equilibrium Flight 21 

4. Airplane in Disturbed Flight 21 

5. Block Diagram of Plant and Controller 34 

in Linear Tracking Problem 

6. General Desired Altitude Trajectory 46 

7. Desired Altitude Trajectory 53 

8. Desired Revised Altitude Trajectory 53 

9. Desired Velocity, Angle of Attack, Pitch 54 

Angle, and Pitch Rate Trajectory 

10. The Solution of the Riccati Equation 67 

for Case I 

11. The Auxiliary Function s^ for Case I 68 

12. The Optimal Control for Case IA 69 

13. The Optimal and Desired Altitude 70 

Trajectories for Case IA 

14. The Optimal 0L, 0, and 0 Histories 71 

and n, h Traces for Case IA 

15. The Optimal Control for Case IB 72 

16. The Optimal and Desired Altitude 73 

Trajectories for Case IB 

17. The Optimal a, 0, and 0 Histories 74 

and n, h Traces for Case IB 



5 



Figure Page 

18. The Optimal Control for Case IC 75 

19. The Optimal and Desired Altitude 76 

Trajectories for Case IC 

20. The Optimal O', 0, and 0 Histories 77 

and n, h Traces for Case IC 

21 The Solution of the Riccati Equation 78 

for Case II 

22. The Auxiliary Function j3 for Case II 79 

23. The Optimal Control for Case HA 80 

24. The Optimal and Desired Altitude 81 

Trajectories for Case IIA 

25. The Optimal v, 0i, 0 , and 0 Histories 82 

and n, h Traces for Case IIA 

26. The Optimal Control for Case IIB 83 

27. The Optimal and Desired Altitude 84 

Trajectories for Case IIB 

28. The Optimal v, Ol, 0 , and 0 Histories 85 

and n, h Traces for Case IIB 

29. The Optimal Control for Case lie 86 

30. The Optimal and Desired Altitude 87 

Trajectories for Case IIC 

31. The Optimal v. O', 0, and 0 Histories 88 

and n, h Traces for Case IIC 



6 



ACKNOWLEDGEMENT 



The author is indebted to Professor D. E. Kirk for 
originally suggesting the problem and for providing guid- 
ance, assistance, and encouragement throughout the inves- 
tigation. Appreciation is also extended to Professor 
E. R. Rang for his valuable assistance in aerodynamic 
analysis and to LCDR J. R. Wilson, Jr., who provided the 
majority of the F-4J airplane data. 



7 



I. 



INTRODUCTION 



A most critical portion of flight, even in an ideal 
operating environment, is the landing of an airplane at a 
desired landing site. Inclement weather at the landing 
site, with the corresponding reduction in visibility, fur- 
ther complicates this portion of the flight operation to 
the point that landing may become an impossibility due to 
visibility restrictions. The increased demand, both com- 
mercially and militarily, for flight operations under all- 
weather conditions has prompted a great deal of research 
to develop an "Automatic All-Weather Landing System" capa- 
ble of automatically controlling an airplane during the 
landing portion of a flight. The problems encountered in 
such investigations, although not considered insurmount- 
able, are vast and complicated mainly because the flight 
dynamics of the airplane and the surrounding environment 
are themselves complicated. This investigation makes some 
simplifying assumptions to reduce this complexity and de- 
rives the necessary controls for such an all-weather land- 
ing system, using optimal control techniques. 

The problem and the simplifying assumptions are de- 
scribed in Chapter II. In Chapter III a general model for 
an airplane in the landing mode is developed, and the spe- 
cific model used in this investigation is presented. 
Chapter IV presents the optimal control theory employed 
and describes the method used to modify the original state 
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variables to make the theory applicable. The problem 
specifications are introduced in Chapter V, the desired 
trajectory is formulated, and the significance of the per- 
formance measure in obtaining realistic results is dis- 
cussed. In Chapter VI the procedures used in the investi- 
gation are explained, and the cases investigated are pre- 
sented. The results are discussed in Chapter VII, and the 
conclusions and recommendations are presented in Chapter 
VIII. 
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II. PROBLEM FORMULATION 



The basic descriptions and definitions of the landing 
portion of flight are presented in this chapter, along with 
the assumptions made to reduce the complexity of the pro- 
blem. 

A. STANDARD WEATHER CRITERIA 

The International Civil Aviation Organization (ICAO) 
has adopted the following set of minimum weather conditions 
for the automatic landing of an airplane, (Ref. 1) . All 
distance specifications are given in meters, followed in 
parenthesis by approximate values in feet. 

Category I. Operation down to minima of 60 meters 
(200 ft.) decision height (altitude) and Runway Visibility 
Reading (RVR) of 800 meters (2600 ft.). 

Category II. Operation down to minima of 30 meters 
(100 ft.) decision height and a RVR of 400 meters (1200 ft.). 

Category III-A. Operation to and along the surface of 
the runway, with external visual reference during the final 
phase of the landing to a RVR minimum of 200 meters (700 
ft.) . 

Category III-B. Operation to and along the surface of 
the runway and taxiways, with visibility sufficient only 
for visual taxiing comparable to an RVR of about 50 meters 
(150 ft.) . 
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Category III-C. Operation to and along the surface of 
the runway and taxiway without external visual reference. 

The Federal Aviation Agency (FAA) has established the 
criteria for certification of any automatic all-weather 
landing system based on compliance with these categories of 
automatic landing weather minima. Presently no automatic 
systems that comply with any Category III minima have been 
certified. However, experimental tests have been success- 
fully conducted (Category III-B) with a system in a C-141 
airplane (Ref. 2) . Several automatic control systems have 
been certified for Category II operation, provided that the 
landing site has the prerequisite ground equipment. The 
vast majority of automatic landings are restricted to 
Category I conditions. 

B. GROUND EQUIPMENT AT AN ALL-WEATHER LANDING SITE 

At present, only one special piece of ground equipment 
an Instrument Landing Systems (ILS) , is required for land- 
ing at an all-weather landing facility. The ILS provides a 
radio beam to establish a nominal glide path to the land- 
ing site and several locator beacons, placed along the land 
ing track to provide distance checks. One locator beacon 
is positioned to provide an indication of the nominal deci- 
sion altitude point for the facility. Both azimuth and 
elevation information are provided by the ILS. Category I 
and Category II automatic landing operations are presently 
using this ground equipment with the latter, of necessity. 
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requiring more sophisticated ILS equipment to provide the 
increased accuracy needed in the radio beam. Figure 1 
shows present requirements for the beam accuracy of ILS 
Category-II ground equipment. Several major U.S. air- 
fields have this equipment installed and operational. 
Because of beam distortions due to the electromagnetic dis- 
turbances, antenna installation environments, etc., the ILS 
equipment is not accurate below an altitude of about 60 
feet. Therefore, any automatic control of the airplane 
beyond Category II minima will require either additional 
ground equipment or an automatic system in the airplane it- 
self to effect the final phases of the landing. 




Figure 1 

Category II ILS Beam Accuracy Requirements 
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c. 



PHASES OF ALL-WEATHER LANDING 



The automatic landing of an airplane under all-weather 
conditions can be divided into five phases: 

1. The Approach-Intercept Phase 

The Approach-Intercept Phase is defined as that 
segment of the landing during which the configuration of 
the airplane is physically changed in anticipation of the 
landing and the airplane is directed toward the ILS beam 
which serves the landing site. 

2 . The Capture-Track Phase 

The Capture-Track Phase is defined as that segment of 
the landing during which the airplane is directed toward 
the landing point by an automatic flight control system 
which uses the ILS beam as its reference. 

3 . The Flare Phase 

The Flare Phase is defined as that segment of the 
landing during which the airplane attitude is changed in 
preparation for the touchdown. 

4 . The Land Phase 

The Land Phase is defined as the actual touchdown 
of the airplane on the landing surface. 

5 . The Roll-Out Phase 

The Roll-Out Phase is defined as the climax of 
the landing during which the airplane is slowed to a stop 
or, equivalently, to a safe taxi speed. 
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Figure 2 depicts the phases of landing in both azimuth 
and pitch. Included, for clarity, is the decision altitude 
(height) , which is defined as the altitude at which a deci- 
sion must be made as to whether it is feasible to continue 
the landing. As an example, the decision altitude for 
Category II landings is 100 feet above ground level. 

Automatic control of the airplane during the first 
three phases of the landing is presently standard operating 
procedure for both commercial and military flight opera- 
tions. However, lack of certified automatic flight control 
systems to accomplish the remaining phases of the landing 
dictates that the decision to continue the landing beyond 
the decision altitude be based on the premise that the 
pilot has visual contact with the landing site at the de- 
cision altitude, so that the airplane can be controlled 
manually, if necessary. 

D. ASSUMPTIONS MADE FOR THE LANDING PROBLEM IN THIS STUDY 

As stated in the opening remarks of this section, some 
simplifying assumptions are made in this study to reduce the 
complexity of the automatic landing problem. These assump- 
tions are listed below. 

Assumption 1. Only that portion of the landing from 
the decision altitude to the land phase will be considered. 
It is assumed that the airplane is automatically or manually 
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FIGURE 2 

PHASES OF 
AUTOMATIC LANDING 
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controlled, using the ILS as a reference, to the locator 
beacon which defines the nominal decision altitude point -- 
taken to be 100 feet above ground level. The roll-out 
phase is also assumed to be automatically or manually con- 
trolled. 

Assumption 2. The airplane is physically located in 
space within the prescribed Category-II ILS "window" at the 
nominal decision altitude point and is in Equilibrium 
Flight (see definition in Chapter III) at that time. If 
these conditions are not met, it is assumed that the land- 
ing will be discontinued (airplane will be waved off) . 

Assumption 3. Wind effects will not be considered. 
During landing, the airplane is subjected to both steady- 
state and gusty winds, the latter being of primary impor- 
tance, since the gusts are random in nature,, Steady-state 
winds could be considered since their effect can be elimi- 
nated by a steady-state change in the airplane heading, 
whereas the random wind gusts could require a statistical 
handling of the problem. By neglecting both effects, the 
airplane velocity becomes equal to the ground velocity and 
the problem can be formulated in terms of time-to-go-to 
landing. 

Assumption 4. Only airplane motion in the vertical 
plane will be considered. Lateral motion of the airplane 
in the final phases of landing is primarily necessary to 
point the airplane in the direction of the runway just prior 
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to the actual touchdown. The airplane is normally "off" 
heading during the landing to counter steady-state cross- 
winds. Since the wind effects have been neglected, the 
lateral motion can be neglected. 
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III. MODELING THE PLANT 



The first step in any control problem is the formula- 
tion of a realistic mathematical model to represent the 
dynamics of the plant to be controlled — in this case, the 
airplane. In this chapter the formualtion of a general 
model applicable to any airplane is discussed, and the spe- 
cific model used in this investigation is presented. 

A. THE GENERAL MATHEMATICAL MODEL 

Aerodynamists have developed two sets of linear equa- 
tions which describe the dynamics of an airplane. These 
are referred to as the longitudinal or symmetric and the 
lateral or asymmetric equations of motion. Since, as pre- 
viously stated, this problem considers only the longitudi- 
nal motions of the airplane, further reference is limited 
to the longitudinal equations of motion only. Derivation 
of these equations is beyond the scope of this study, but 
complete and detailed derivations can be found in the li- 
terature. (Ref. 3, 4, and 5) 

Before presenting the equations of motion, it is sig- 
nificant, for clarification and reference purposes, to pre- 
sent several aerodynamic definitions and assumptions used 
in obtaining these equations. 

1 . Aerodynamic Definitions 

a. Equilibrium flight is defined as unaccelerated 
flight. Flight is along a straight path during which the 
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linear velocity vector measured relative to fixed space 
is invariant and angular velocity is zero. 

b. Steady flight is defined as flight during 
which the linear velocity vector is invariant and angular 
velocity is constant. In this context, equilibrium flight 
is also steady flight. 

c. Airplane coordinates, axes, and angles are 
defined in Figure 3, which shows the airplane in equili- 
brium flight. 

v q = Equilibrium flight linear velocity along x axis 
= Equilibrium flight angle of attack 
= Equilibrium flight pitch angle 
y Q = Equilibrium flight glide angle 

d. Instantaneous components are defined as the 
summation of equilibrium flight components and correspond- 
ing perturbations caused by a disturbed flight condition. 
Figure 4 is a representation of the airplane in a disturbed 
flight condition and includes the perturbations of the two 
control components. 

V = Instantaneous linear velocity along x axis 

W = Instantaneous linear velocity along z axis 

Q = Instantaneous angular velocity along y axis 
A = instantaneous angle of attack 

® = Instantaneous pitch angle 
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r = Instantaneous glide angle 

6^ = Elevator deflection perturbation 

6^ = Thrust perturbation 

2 . Assumptions Used in Obtaining Eguations of Motion 

a. The airframe is a rigid body (no aeroelastic 
deflection of the airframe) . 

b. The earth is fixed in space, and the earth's 
atmosphere is fixed with respect to the earth. 

c. The mass of the airplane remains constant 
during any particular dynamic analysis. 

d. The x-z plane (vertical plane) is a plane of 

symmetry. 

e. Disturbances from steady flight are small 
enough so that the products and squares of changes in ve- 
locities are negligible in comparison to the changes them- 
selves, and the disturbance angles are small enough so that 
the sines are equal to the angles in radians and the cosines 
are equal to one. Products of these angles are also approx- 
imately zero. Because these disturbances are small, the 
change in air density encountered by the airplane during 
any disturbance is considered to be zero. 

f. During steady flight conditions, the airplane 
is assumed to be flying with wings level and all components 
of velocity zero except for the linear velocity along the x 
axis. (Equilibrium flight) 
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g. The flow is quasi-steady. (The air flow pat- 
tern around the airframe instantaneously changes in a 
steady flow pattern as the airplane changes its orienta- 
tion with respect to its flight path.) 

3 . The Longitudinal Equations of Motion 

v(t) = X v v(t) + X^q(t) - g9(t) + x^ictt) 

+ X^W(t) + X 6 6 (t) + X 6 6 (t) 
e T 

05 (t) = Z v v(t) + V Q q(t) + Z q(t) + Z w «(t) 

(3.1) 

+ Zj05(t) + Z 6 6 e (t) + Z 6 6 T (t) 
e T 

q(t) = M v v(t) + M q(t) + J^,tO( t) + M^d>(t) 

+ m 6 6 (t) + m 6 6 (t) 

e e °T i 



where the subscripted X, Z, and M are the Dimensional Sta- 
bility Derivatives of the airplane which are parameters for 
a specific airplane in a specific flight regime (e.g., 
landing regime, subsonic regime, etc.). 

By substituting the relationships 



o)(t) = v Q a(t) 
q(t) = 0(t) 



(3.2) 
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the equations of motion can be rewritten in the more fa- 
miliar terms of v, CL, and 0 which follow: 



v(t) = x v v(t) - g0(t) + x q 0(t) + x a a(t) 



+ X^a(t) + X 6 6 (t) + X, 6 (t) 

e T 



z z z 

&(t) = — ^ v(t) + 0(t) + ~ 0(t) + —■ a(t) 
o o o 



+ — &(t) + — 

V V 

o o 



6 + — - 6 
e v T 



( 3 . 3 ) 



• • • 

0(t) = M v v(t) + M q 0(t) + M a a(t) + M & a(t) 

+ Ms 6 (t) + M 6 6 (t) 

In a rigorous mathematical sense, these longitudinal 
equations of motion, which deal with perturbations from 
equilibrium flight, are applicable only to infinitesimal 
disturbances; however, aerodynamic experience has shown 
that quite accurate results can be obtained by applying 
these equations to disturbances of finite, non-zero 
magnitude. 

Since this problem deals with the landing approach, 
two additional parameters must be considered. One, the 
altitude of the airplane, is of primary importance, and an 
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equation relating actual airplane altitude to the flight 
dynamics is required. Such a relationship can be approx- 
imated by 

h(t) = v Q r(t) (3.4) 

where h is defined as the instantaneous vertical distance 
above the ground of the aircraft wheels. By definition 
(see Figure 4) 

r(t) = y + y(t) 

° (3.5) 

y(t) = e ( t) - a(t) 

and therefore, 

h(t) = v Q y Q + v Q (9(t) - a( t) ) (3.6) 

The second parameter is ground effect — defined as the 
effect on the airplane dynamics as the airplane nears the 
vicinity of the stationary plane (the ground). This ef- 
fect on an airplane is non-linear and is dependent on the 
size and shape of the specific airplane. 

The formulation of a general mathematical model to 
describe the motions of any airplane in the landing ap- 
proach is now complete. All that remains is application 
to a specific airplane. 

B. THE SPECIFIC MATHEMATICAL MODEL 

The main problem associated with the use of the gener- 
al mathematical model is determination of the numerical 
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values for the Dimensional Stability Derivatives for the 
specific airplane. Because of the vast amount of experi- 
mental data available on the airplane, both from wind tun- 
nel studies and actual flight testing, the McDonald 
Douglas F-4J Phantom II jet fighter was selected. Appendix 
A contains representative data for this airplane in the 
landing configuration. 

As the F-4J Dimensional Stability Derivatives, X , Xi, , 

J q a 

Z q, and M 6 are zero an< ^ ground effects are neglected, 

the longitudinal equations of motion and the height equa- 
tion can be simplified as follows: 




« T 

« e (t) + (t) 

o 



o 



o 



o 



(3.7) 



6(t) = M v v(t) + M^(fc) + M 4 - M 0(t) + M § 6 e (t) 



h(t) = V Q (0(t) - Q!(t)) + V Q y o 



Defining 



x 





6 



) 



(3.8) 
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and 



u 




(3.9) 



the equations are manipulated into the state variable form 



x(t) = Ax( fc ) + Bu(t) + C 



(3.10) 



with the result 




(3.11) 
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where 



a, , = X 
11 v 



a c „ = - v 
52 o 



a 12 x a 



*53 



13 



-g 

Z 



v 



21 v 



12 



21 



= X, 



Z 6< 



a, 



22 v 



M*Z 
a v 



41 



42 



M + 
v v 



M + ^ 

M a + ~ v 



22 



41 



42 



M « Z 6 
m * +-1H 

e o 






v 



a. ., = M + M* 
43 q a 



c _ i = v y 
51 o o 



Appendix B contains the numerical values for these elements. 

Heretofore in the formulation, it was assumed that the 
two perturbation controls, elevator deflection and thrust, 
are instantaneously available to control the airplane. In 
reality, both controls have associated dynamics, which are 
aerodynamically referred to as elevator actuator lag and 
lag between the thrust command and thrust. These effects 
can be approximated by use of the following differential 
equations : 



6 (t) = (6 (t) - 6 (t)) 

6 c 

e (3.12) 

« T (t) = (« (t) - 6 (tn 

6 t 
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where the subscripted T's are the appropriate component 

time constants and 6 and 6 m are the elevator deflection 

e T 

c c 

and thrust commands respectively. In this investigation 
it will be assumed that the controls are instantaneously 
available; i.e., the dynamics will be neglected. 

This completes the mathematical modeling of the F-4J 
airplane in the landing configuration. The model was test- 
ed on an IBM System/360 digital computer, using a Fourth- 
Order Runge-Kutta routine to solve the differential equa- 
tions. The responses resulting from various control in- 
puts were consistent with actual airplane responses. 
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IV. OPTIMAL CONTROL THEORY 



This chapter will not attempt to present the entire 
field of optimal control theory, but rather will assume 
that the reader is familiar with it; hence, only those 
techniques pertinent to this investigation will be dis- 
cussed. Some definitions will first be presented, follow- 
ed by mathematical techniques employed to obtain the op- 
timal control. 

A. OPTIMAL CONTROL DEFINITIONS 

1 . Control History 

The history of control input values during the 

time interval [t , t_] is denoted by u. 

or — 

2 . State Trajectory 

The history of state values during the time in- 
terval [t , t,.] is denoted by x. 

o f J — 

3 . Admissible Control 

U denotes the set of all control histories which 

satisfy the physical control constraints during the time 

interval [t , t ,.] . 

o f 

4 . Admissible Trajectory 

X denotes the set of all state trajectories which 

satisfy the physical state variable constraints during the 

time interval [t , t,.] . 

o r 
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5. The Performance Measure 



J denotes a scalar measure of the performance of 
a system when a control history is applied. 

6. Optimal Control 

A control which causes the system to follow an 
admissible trajectory which minimizes the performance mea- 
sure is called an optimal control and is denoted by u*. 

7 . Optimal Trajectory 

The admissible trajectory which results when an 
optimal control is applied to the system is called an op- 
timal trajectory and is denoted by x*. 

8 . Tracking Problem 

A problem wherein the intent is to maintain the 
state trajectory as close as possible to a desired trajec- 
tory — denoted by _r — in the interval [t^t^] is called 
a tracking problem. 

B. THE LINEAR TRACKING PROBLEM 

Since the automatic landing problem as formulated in 
Chapter II is directed toward controlling the airplane in 
a desired manner during the final phases of a landing, and 
the mathematical model of the airplane as presented in 
Chapter III is linear, the problem can be considered as a 
linear tracking problem. 
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It has been shown (Ref. 6 and 7) , where complete deri- 
vations can be found, that, given a set of state equations 



x(t) = A(t) x(t) + B(t)u(t) 



(4.1) 



where 

x(t) is the state vector of dimension n 
u(t) is the control vector of dimension m 
A ( t) is a nxn matrix 
B ( t) is a nxm matrix 
and the performance measure 




2 

H 



p r _ 



+ 2 



x(t) - r (t) || + ||u(t) || Jdt 

Q(t) R ( t) 



(4.2) 



where 



l|x( T ) 
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the final time t^ is fixed, x(tf) f ree » u(t) is uncon- 
strained, H and Q(t) are real symmetric positive semi- 
definite matrices, and R(t) is a real symmetric positive 
definite matrix, that the optimal control exists and is 
unique. The optimal control is given by 

u*(t) = F ( t) x*(t) + g(t) (4.3) 



32 



where F(t) is the mxn matrix of feedback gains, and g(t) 
is the mxl command signal vector which is dependent on 
the system parameters and the desired state trajectory. 
Figure 5 is a block diagram of the plant and optimal con- 
troller. 

F(t) and g(t) are given by 



F ( t) = - R ( t) 1 B T (t)K(t) (4.4) 

g(t) = - R ( t) -1 B T (t) s(t) (4.5) 



where K(t) is the solution of the Riccati-type matrix dif- 
ferential equation 



K(t) = - K(t)A(t) - A T ( t) K ( t) - Q ( t) 
+ K(t)B(t) R -1 ( t) B T ( t) K ( t) 



(4.6) 



with boundary conditions 



K(t f ) = H (4.7) 

and s^(t) is the solution of the linear vector differential 
equation 

i(t) = - [A T (t)-K(t)B(t)R _1 (t)B T (t) ]s(t)+Q(t)r(t) 

(4.8) 

with boundary conditions 

s.(t f ) = - Hr(t f ) (4.9) 
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PLANT 




FIGURE 5 

BLOCK DIAGRAM OF PLANT AND 

CONTROLLER 

IN A LINEAR TRACKING PROBLEM 
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Since the boundary conditions at the final time t f are 
known, these two matrix differential equations must be in- 
tegrated backwards in time from t to t . 

f o 

It was assumed that all states of the system can be 
measured. 

C. MODIFICATIONS TO SPECIFY STATE PLANT 

The previous discussion indicated that the optimal con- 
trol in a linear tracking problem exists and can be obtain- 
ed, assuming that the plant is of the form 

x(t) = A ( t) x ( t) + B ( t) u ( t) (4.10) 

Since the plant model, as presented in Chapter III, is in 
the form 

x(t) + Ax(t) + Bu(t) + c_ (4.11) 

some modifications are necessary in order to adapt the 
specific problem to the theory presented. 

This modification can be accomplished by observing 
that the c_ vector contains only one non-zero element, i.e., 
the differential equation describing the altitude of the 
airplane 

h ( t) = v Q y o + v Q (0(t) - a(t)) (4.12) 

contains the constant term v y 

o o 
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By defining 



h(t) = h(t) - h e (t) (4.13) 

where 

h g (t)= h (0) + v o y Q t (4.14) 

is the nominal equilibrium altitude of the airplane during 



the interval [t , t-] . It follows that 
o f 

K {t) = v o 7 o (4 * 15) 

h(t) = fi(t) - h e (t) (4.16) 

h(t) = v o r o + v Q (0(t) - a( t) ) - v Q y o (4.17) 

h ( t) = v o (0(t) - a(t)) (4.18) 



So h(t), as defined, is the perturbation of the altitude of 
the airplane about the nominal equilibrium altitude. 

With this modification, a new set of states, hereafter 
referred to as the revised states (to distinguish them from 
the actual states) , can be defined as 



x 





(4.19 
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and the revised state equations become 



x(t) = Ax ( t) + Bu ( t) (4.20) 

where the elements of A and B remain identical to those pre- 
viously formulated in Chapter III, equation (3.11). This 
revised plant is now in the proper form for application of 
the optimal control theory. The solution then becomes 

u*(t) = F(t)x*(t) + g(t) (4.21) 

from which the behavior of the actual states can be deter- 
mined by substituting the definition of x(t). This substi- 
tution produces the optimal control 




(4.22) 



By this simple mathematical manipulation, the problem can 
be solved by applying optimal control theory to the revised 
state equations and then expressing the results in terms of 
the actual states. 
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V. SPECIFICATIONS, DESIRED TRAJECTORIES 



PERFORMANCE MEASURE 

In this chapter the applicable specifications are in- 
troduced and a set of desired state trajectories for the 
problem are derived. Many of the results presented are 
based on the author's experience as a pilot operating the 
F-4 and other jet airplanes. The final section of this 
chapter discusses the significance of the performance meas- 
ure in obtaining realistic results. 

A. PROBLEM SPECIFICATIONS 

A successful automatic landing requires that certain 
conditions relative to the airplane and its environment be 
satisfied. These conditions, which assume the role of spe- 
cifications for the problem, are formulated in terms of per- 
formance requirements and constraints on the system states 
and controls, together with specifications regarding related 
parameters. The following specifications are of primary 
importance to the problem under consideration. 

1. Velocity 

The instantaneous velocity of the airplane during 
the landing portion of flight must remain above 1.1 
or V(t) s 214.5 ft/sec. The upper velocity limit depends 
on the structural limit with landing gear and flaps extend- 
ed; flight experience indicates that a reasonable upper 
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limit is 1.1 v q or V(t) ^ 245 ft/sec. In terms of the 
specified state, v(t), which is the perturbation velocity 
about equilibrium flight, these limitations are 

-8.5 ft/sec s v(t) ^ 22 ft/sec , tc[t Q ,t f ] . 

2. Angle of Attack 

The angle of attack must remain below 0.9 & . , n 

stall 

or X(t) ^ 0.44 radians. No stringent lower limit exists 
but, as in the velocity case, experience indicates that a 
reasonable minimum value is approximately 0.750^ or 
K(t) * 0.25 radians. As a result, the perturbation angle 
of attack limits were established as 

I a(t) | £ 0.11 rad , t*[t ,t f ] . 

3 . Pitch Angle 

During the landing portion of flight, the instan- 
taneous pitch angle is closely related to the angle of 
attack but does not have stringent limitations on its ex- 
cursions. However, at the actual touchdown point, limita- 
tions do exist primarily to prevent the airplane from 
either landing nose wheel first (low pitch angle) or tail 
first (high pitch angle) . For the F-4J airplane these lim- 
itations, in terms of the perturbation pitch angle from 
equilibrium flight conditions, were established as 

-0.20 rad ^ 0(t) =£ 0.25 rad , te[t , t^3 . 
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4 . Pitch-Angle Rate 



Flight experience has indicated that this para- 
meter should be held to a minimum for pilot comfort. (See 
discussion under Related Specifications, below.) Realistic 
limits for the pitch angle rate were established as 

I 0(t)| ^ 0.08 rad/sec , te[t o , t f ] . 

3 . Altitude 

In the problem formulation. Chapter II, it was 
assumed that the altitude of the airplane at the initial 
time was within the prescribed Category II ILS "window" or 

88 ft ^ h ( t ) £ 112 ft 

v o' 

In terms of perturbations about equilibrium flight condi- 
tions, this limitation can be expressed as 

| h(t Q )| ^ 12 ft 

As the airplane approaches the actual touchdown point, how- 
ever, altitude excursions about desired conditions should 
satisfy more stringent limits. As a result, the following 
limitations were established for altitude perturbations 
about a desired altitude trajectory, in the intervals in- 
dicated. (See Section B, this chapter, for desired trajec- 
tory. ) 

In the interval from the decision altitude to the 
flare initiation altitude 

I h(t) - h d (t)| ^ 12ft , t€ C t Q , t-L 3 , 
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and in the interval from the flare initiation altitude to 



actual touchdown 

I h(t) - h d (t)l £ 5ft, tc[t 1 ,t f ] , 

where h d (t) is the desired revised altitude state trajec- 
tory, and t^ is the time at which the flare phase begins. 

6 . Elevator Deflection 

Physical limits for elevator travel exist, and for 
the F-4J airplane in equilibrium flight in a landing ap- 
proach, were established as 

-0.26 rad ^ 6 e (t) s 0.22 rad, tc[t 0 ,t^] . 

7 . Thrust 

Physical limits on thrust available exist for the 
F-4J airplane and are dependent on the selection of power; 
i.e., MILITARY engine operation or AFTER BURNER (A/B) en- 
gine operation. Flight experience has shown that realistic 
limits for thrust perturbations about equilibrium flight 
conditions in the landing portion of flight can be estab- 
lished as 

| 6 T (t)| ^ 3000 lbs, te[t o ,t f ] . 

This range is well within the thrust-available range of 
the F-4J airplane without A/B selection. 
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8 . Related Specifications 



Several related performance requirements and con- 
straints are applicable to the problem and are presented 
in the following listing: 

a. The actual touchdown point must be within 150 
feet of the desired touchdown point on the landing surface. 
Since the problem was formulated in the time domain, this 
restriction was established by requiring that the actual 
touchdown occur within ± 0.65 seconds of the desired touch- 
down time . 

b. At the actual touchdown time the airplane must 
have no tendency to float, or, in other words, the airplane 
should have a positive rate of descent. In addition, the 
rate of descent at touchdown must be within the structural 
sink rate limits for the airplane (Appendix A) . Experience 
indicates that realistic rate of descent limits for the F-4J 
at the actual touchdown time are 

-9ft/sec ^ h(t f ) ^ -3ft/sec 

c. The normal acceleration'*' of the airplane — 
defined as a perturbation from the gravitational accelera- 
tion of equilibrium flight (l.Og) — can be approximated by 

n(t) = — p(6 ( t) - a(t) ) 



1 

The airplane acceleration perpendicular to the hori- 
zontal reference plane. 
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and has physical limits based on the structural strength of 
the airplane. These limits are normally not attained in 
the landing portion of flight, due to the low velocity and 
control effectiveness in this regime. However, the func- 
tion n(t) gives an indication of the smoothness of the land- 
ing, and its effect can actually be felt by the pilot during 
the landing: therefore, to maintain the normal accelera- 

tion as closely as possible to the equilibrium flight condi- 
tion of l.Og the following limitation was established 

| n( t) | ^ 0.2g, tc[t Q ,t f ] . 

B. DESIRED TRAJECTORY 

The problem has been formulated as a tracking problem; 
hence the desired state trajectory must be determined. This 
section presents a general desired trajectory appropriate 
for any airplane in the landing portion of flight, and in- 
dicates how these trajectories are applied to the F-4J 
landing problem. 

1 . General Desired Trajectory 

Since the altitude of the airplane is of great 
importance during the landing, it is an obvious starting 
point in the formulation of a set of desired state trajec- 
tories. In Chapter II it was stated that this investiga- 
tion would consider only that portion of the landing from 
the decision altitude (100 feet) to the touchdown point. 
Referring to Figure 2, this portion of the landing includes 
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a section of the capture-track phase, the flare phase, and 
the land phase. In Chapter III it was shown that the pro- 
blem could be formulated in the time domain. These two 
considerations lead to a general altitude trajectory (spe- 
cified in the time domain) consisting of a constant descent 

from the decision altitude at t = t = 0 to a selected 

o 

flare point at t = t^, followed by a flare to touchdown at 
t = t f . 

In the optimal control theory presented in Chapter IV, 
one of the necessary requirements was that the final time 
(tf) be fixed. To alleviate the obvious disadvantage that 
this requirement imposes, it seemed appropriate to extend 
the final time beyond the actual touchdown time to ensure 
that a landing occurs prior to t f . By so doing, an imag- 
inary flare plane below the actual landing plane was 
established. This implies that the airplane nominally will 
land before the final time specified. Figure 6 depicts 
this general concept in graphical form, where the capture- 
track phase is the interval [o,tf), the actual flare phase 
is the interval [tf,t 2 )> and the actual land phase is t£. 
The imaginary flare phase is the interval [t^,tf) and the 
imaginary land phase is at the time t^. 

The rate of altitude change (rate of descent), h(t), 
must be considered in conjunction with the general altitude 
trajectory, even though it is not a state variable. Ideal- 
ly, during the capture-track phase [o,t^), the rate of 
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descent will be a constant, while during the flare phase, 
it will steadily decrease so that at the actual touchdown 
time (^ 2 ) , the rate of descent will be within the estab- 
lished limits. Utilizing this consideration, the set of 
equations 



h d (t) = h e (t) + h (0) + v Q y o t (5.1) 

h a (t) = h e (t) = v o y o (5.2) 

which describe the desired altitude trajectory, h d (t), and 
rate of descent, h d (t), for the interval [0,t^) were for- 
mulated. 

For the interval [t. ,t J _3 , the desired altitude and 

1 r 

rate of descent were defined as 



L 2 r L i (t_t i ) -\ 

h d (t) = h (0) + v o y Q t + -lj-L^t-t^ 

(5.3) 

f L (t-t ) . 

h d( fc ) = Vo + L 2 ( e -0 ( 5 - 4 > 



Since the revised altitude state was defined in Chapter IV 
to be 

h(t) = h(t) - h e (t) (5.5) 

equations (5.1) and (5.2) became 

h d (t) = h d (t) = 0, te[0,t) , (5.6) 
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Figure 6 

General Desired Altitude Trajectory 



and 



h d (t) 



0 



i -1) - L 



2 (t t x ) 



(5.7) 



r s L-, ( t- 1. ) s. 

h d (t) = L 2 (e x -1 ) , t€[t l ,t f ] . (5.8) 

Equations (5.6) and (5.7) represent the desired revised 
altitude trajectory in the interval [o,t^] in terms of per 
turbations about the equilibrium altitude. 
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These general expressions for the desired revised rate 
of descent, h^(t), provide insight into the form of the 
desired angle of attack trajectory, 0^, and pitch angle 
trajectory, © d , since from equation (4.18) 

h ( t) = v o (9(t) - a(t) ) . (5.9) 

Therefore, 

6 d (fc) " a d (fc) = °' t^Co,^) (5.10) 

and 

L~ s L. (t— t. ) v 

e d (t) - <* d (t) = ^ (e 1 1 -1) , t*[t ,t f ], (5.11) 

O 

By definition (Chapter II) , 0 1 and 0 are zero when the air- 
plane is in equilibrium flight, hence equation (5.10) is 
consistent. The problem remains to find equations for a 
and 0 in the interval [t^,t^] which are consistent with 
the requirements of equation (5.11). Before proceeding, 
recall that the desired pitch-angle-rate trajectory, 0 
must also be considered. An acceptable expression for the 
desired pitch-rate trajectory is a linear function during 
the flare phase interval [t^t^] given by 

0 d (t) = L^t-t^ (5.12) 
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Integrating equation (5.12) and applying the desired bound- 



ary condition at t = t^, i.e., = 0 



9 a (t> = -r (t - fc i )2 - 



(5.13) 



Substituting equation (5.13) into (5.11) yields 



m ^3 , _ . 2 _ ^2 / L l (t t l ) _ n 

tt d^ 2 ^ ~ fc l^ v v e 1 y 



(5.14) 



as the expression for the desired angle of attack in the 
interval [t^,t^]. 

The final state component trajectory which must be 
considered is the velocity trajectory. Experience from 
both an operational and a safety standpoint has establish- 
ed that the airspeed during a landing should remain es- 
sentially constant. This implies that v^, the desired 
velocity trajectory, should be zero during the entire in- 
terval of interest, hence 



v d (t) = 0, tc[0,t f ] . (5.15) 

2 . Specific Desired Trajectories 

In the preceding section a general set of desired 
state trajectories — applicable to any airplane — were 
formulated. What remains is application to the F-4J air- 
plane so that the unknown constant terms can be evaluated. 
Since, in the interval [0,t d ), the desired state trajec- 
tory is zero, the following discussion will be concerned 
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with the interval [t^,t^] where the equations of interest 
are (5.3), (5.4), (5.7), (5.8), (5.12), (5.13) and (5.14) 

— hereafter referred to as the reference equations. 

The first consideration, however, was to determine a 
realistic time frame for the problem. The equilibrium air- 
speed, v , of the F-4J, as given in Appendix A, is 223 
ft/sec. Assuming that the airplane is at the nominal de- 
cision altitude (100 feet) at t = 0 and that the airplane 
remained in equilibrium flight on a negative three-degree 
glide slope, the actual time of touchdown from equation 
(4.14) is approximately 8.57 seconds. This assumes that 
there is no flare phase. Although the standard operating 
procedure in the U.S. Navy is to operate the F-4J airplane 
with no flare, this investigation will include a flare 
which commences at t^ = 6 seconds; at this time the air- 
plane is at a nominal altitude of approximately 30 feet. 

It follows then that a reasonable time frame for the prob- 
lem is 10 seconds, and that the interval [t^,t^] is 4 sec- 
onds in duration. 

Next, values for the unknown constants in the refer- 
ence equations must be obtained. Boundary conditions for 
the reference equations at t = = 6 seconds are known, 

and final-time boundary conditions can be selected to pro- 
vide realistic trajectories. However, recalling that the 
actual touchdown of the airplane occurs prior to the final 
time, and that the specifications refer to this actual 
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touchdown time, another set of boundary conditions at 
t = t 2 must be satisfied. The conditions at t = t 2 are not 
exactly given, but realistic values can be selected to meet 
the necessary specifications. As a result, values for the 
unknown constants L^, L 2 , and of the reference equations 
were found by a trial-and-error method. Values of these 
constants were arbitrarily chosen to satisfy boundary con- 
ditions at t = t^ and then substituted into the reference 
equations to obtain the actual touchdown time, and the 

corresponding values of the parameters at t = t 2 . This 
process was continued until values of L^, L 2 , and were 
found, which resulted in realistic trajectories. The values 

1 j 1 = 0.25 

L 2 = 4.68 

L 3 = 0.016 

which were obtained result in a desired actual touchdown 
time, (t 2 ), seconds. The reference equations which 

apply to the interval [6,10] are then 




(5.16) 




(5.17) 
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(5.18) 



S d (t) = otH (e* 25(t-6) -l) -4.68 (t-6) 

h d (t) = 4.68 (e* 25(fc_6) -l) (5.19) 

e d (t) = 0 . 016 ( t-6 ) (5.20) 

0 d (t) = 0 . 008 ( t-6 ) 2 (5.21) 

Of (t) = 0.008(t-6) 2 - ^e* 25(t_6) -l) , (5.22) 

o 

while in the interval [0,6) the applicable equations are 

h d (t) = h (0) + v Q y o t (5.23) 

h d <« = v o y Q (5.24) 

h d (t) = E d (t) = e d (t) = 6 d (t) = « d (t) = 0 (5.25) 

and 

v, ( t) = 0, tf[0,10] . (5.26) 

d 



Table A summarizes the solutions of equations (5.16) through 
(5.22) at the flare initiation time (t^) , the actual touch- 
down time (t 2 ), and the final time (t f ). 
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Table A 



Solutions to Reference Equations at Flare Initiation 
Time, Actual Touchdown Time, and Final Time 



Solutions at 



Parameter 


Equation No. 


t=6 sec 


t=9.3 sec 


t=10 sec 


h d (ft) 


5.16 


29.947 


0 


-3.286 


• 

h^( ft/sec) 


5.17 


-11.676 


-5.667 


-3.621 


h d (ft) 


5.18 


0 


8.582 


13.469 


• 

h d (ft/sec) 


5.19 


0 


6.009 


8.055 


• 

0 d (rad/sec) 


5.20 


0 


0.053 


0.064 


0 d (rad) 


5.21 


0 


0.087 


0.128 


« d (rad) 


5.22 


0 


0.060 


0.092 



Figure 7 presents a plot of the actual desired altitude 
trajectory with the equilibrium altitude trajectory includ- 
ed for comparison. Figure 8 illustrates the revised desired 
altitude trajectory h^, while Figure 9 depicts the three 
angular components of the desired state trajectory, , 

© d , and 0^ and the velocity trajectory v^. 
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Figure 7 



Desired Altitude Trajectory 




Figure 8 

Desired Revised Altitude Trajectory 
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* 




Figure 9 

Desired Velocity, Angle of Attack, Pitch Angle, and 
Pitch Rate Trajectory 
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C . PERFORMANCE MEASURE 

The mathematical form of the performance measure in- 
troduced in Chapter IV was 

2 

J = 1/2 || x(fc f ) “ r(t ) || 

H 




(5.27) 

One of the advantages of this quadratic performance measure 
is that the elements of the H, Q(t) and R(t) matrices, which 
shall be called the weighting matrices, can be related to 
the design parameters of the system and chosen to satisfy 
the design objectives. If these weighting matrices meet 
the requirements established in Chapter IV, specifically 
H and Q ( t) are positive semi-definite and R(t) is positive 
definite, the optimal control law can be found, and is 
unique for the selected set of weighting matrices. 

Previous investigations dealing with similar airplane 
landing problems (Refs. 8 and 9) have indicated that the 
weighting matrices must be diagonal matrices with positive 
diagonal entries. In other words, the weighting matrices 
must have non-zero weighting values specified for each 
state and control. As a result, in the initial trials it 
was assumed that the weighting matrices are diagonal and 



55 



time invariant over the interval of interest. Since the 
form of the weighting matrices was so selected, each dia- 
gonal element of the respective matrix could be related to 
a specific state variable or control, and therefore, theo- 
retically, could be assigned a proper value to meet estab- 
lished specifications. Unfortunately, this selection of 
weighting factors was not simply accomplished. In fact, 
it became obvious during the investigation that the state 
trajectory was very sensitive to changes in the weighting 
matrices. (See related discussion in Chapter VII.) 

I 

The selection of values for the weighting matrices was 
carried out by a trial-and-error method wherein a set of 
weighting matrices was selected, the optimal control law 
computed, the optimal trajectory calculated, and the re- 
sults compared with established specifications. In effect, 
the final weighting matrix selection was based on obtain- 
ing a realistic optimal trajectory which conforms to the 
specifications . 

Appendix B contains the numerical values used in the 
weighting matrices for this investigation. 
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VI . INVESTIGATION PROCEDURE 



Since the previous discussion has presented the pro- 
cedures used in modeling the plant, developing the desired 
trajectory, and related topics, only the procedure used to 
obtain the optimal control will be presented. The actual 
cases investigated will also be discussed in the last sec- 
tion of this chapter. 

A. PROCEDURE 

The following procedure was used to obtain the op- 
timal control: 

1. A set of representative values were selected for 
the diagonal elements of the weighting matrices (H, Q, and 
R) of the performance measure, equation (5.27). 

2. Equations (4.6) and (4.8) were simultaneously in- 
tegrated from t = tf to t = 0 to obtain values for K(t) 
and s^(t) respectively. 

3. Using the results of step 2, equations (4.4) and 
(4.5) were solved to obtain values for F(t) and g(t) re- 
spec tively . 

4. The state equations (3.11) were integrated from 
t = 0 to t = tf, using the results of step 3 to obtain 
u*(t) of equation (4.22). 

5. The optimal control and optimal trajectory ob- 
tained from step 4 were observed to ascertain compliance 
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with problem specifications. If the results were not ad- 
missible or unrealistic, the values of the elements of the 
weighting matrices were changed and steps 2 through 4 were 
repeated. 

To implement the foregoing procedure, a computer 
program was written in Fortran IV for use in the IBM Sys- 
tem/360 digital computer. System/360 Scientific Subrou- 
tines were used in the program to accomplish the matrix 
algebra. An existing Fourth-Order Runge-Kutta subroutine 
developed at the Naval Postgraduate School was used for 
the required numerical integration. The computer program, 
including the integration subroutine, is presented after 
the appendices. 

B. CASES INVESTIGATED 

Two cases were investigated in the initial study. 

Case I assumes that the velocity of the airplane is con- 
stant and elevator deflection is the only control avail- 
able; Case II considers the plant and controls as pre- 
viously formulated. As a result, two different state mo- 
dels were used. 

1. Case I 

Since it is assumed that the velocity is constant 
and that elevator deflection is the only control available, 
the state model, equation (3.11) reduces to 
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( 6 . 1 ) 




and 



u 




( 6 . 3 ) 
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By the same argument as was presented in Chapter IV, this 
state equation can be manipulated into the revised state 
equation 



x(t) - Ax(t) + Bu(t) 



(6.4) 



where 



x 





a 



0 



0 



) 



(6.5) 



and the elements of A and B remain identical to those of 
equation (6.1). 

2. Case II 

The state model, equation (3.11), and the revised 
state model, equation (4.20), as previously presented, were 
used in Case II. 

In order to evaluate the optimal trajectory, three dif- 
ferent sets of initial conditions were selected. The first 
assumes that the airplane is in an ideal flight condition, 
or, in terms of the revised state model, all initial state 
values are zero. The other two sets were selected to con- 
form to the worst initial flight conditions that could be 
reasonably expected. The following listing discusses 
these initial flight conditions. 
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a. HIGH AND FAST. This set of initial conditions 
assumes that the airplane was at the upper altitude limit 
of the ILS Category II "window" and that the airplane ve- 
locity is above v , the equilibrium airspeed. In Case I, 
the second assumption is not, strictly speaking, applicable 
since it is assumed that velocity was constant. However, 
the condition was simulated by assuming that the initial 
angle of attack, a(t Q ), and pitch angle, 9(t Q ), were below 
their respective equilibrium values. This condition would 
result in actual flight if the airplane velocity were 
above the equilibrium flight velocity. It was further 
assumed that the initial pitch angle rate, 9(t Q ), was zero. 

b. LOW AND SLOW. This set of initial conditions 
assumes that the airplane is at the lower altitude limit of 
the ILS Category II "window" and that the airplane velocity 
is below the equilibrium airspeed. By the same reasoning as 
before, the second assumption was simulated for Case I by 
assuming that the initial angle of attack and pitch angle 
were above their respective equilibrium values. Again, 

the initial pitch angle rate was assumed to be zero. 

Table B presents a summary of the cases investigated 
and includes the initial values used for each case. 
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Table B 



Summary of Cases Investigated 



Case 



Initial Flight 
Condition 



Applicable Initial State 

State Equations Condi tions 



I A IDEAL 

IB HIGH AND FAST 

IC LOW AND SLOW 

II A IDEAL 

II B HIGH AND FAST 

IIC LOW AND SLOW 



(6. 1,6. 4) (0,0,0,0) T 

(6. 1,6. 4) (-0.03, -0.03,0, 112. 0) T 



(6. 1,6. 4) (0.03,0.03,0,88.0) 



(3.11,4.20) 



(0,0,0,0,0) T 



(3.11,4.20) (5. 0,-0. 03, -0 .03, 0,112 JO? 

(3.11,4.20) (-5. 0,0. 03, 0.03,0, 88. 0) T 
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VII. RESULTS 



A general result is first presented, followed by spec- 
ific results, in graphical form, for each of the cases in- 
vestigated. The final section of this chapter discusses a 
problem encountered during the investigation. 

A. GENERAL RESULTS 

The optimal control for an all-weather landing in the 
F-4J airplane was derived (see Specific Results below) 
using the techniques presented. Since some simplifying as- 
sumptions were made to reduce the complexity of the prob- 
lem, the results are not conclusive, but serve to demon- 
strate that the design of an automatic controller for the 
landing of an airplane is feasible by formulating the prob- 
lem as a tracking problem and applying optimal control 
theory. By systematically eliminating some of the assump- 
tions made in this study, a closer approximation to the 
actual airplane landing problem can be accomplished. 

Actual implementation of an automatic controller to 
provide the optimal control derived by these techniques 
would require that a digital computer be placed aboard the 
airplane which may impose an undersirable penalty. If so, 
the optimal control still provides the goal for the de- 
sign of any sub-optimal controller and hence is of great 
value to the control engineer. 
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B. SPECIFIC RESULTS 

Specific results for each of the cases investigated 
are presented in graphical form and are summarized below. 

1. The solution to equation (4.6), the K(t) values, 
are presented for: 

a. Case I in Figure 10 

b. Case II in Figure 21 

2. The solution to equation (4.8), the s^(t) values, 
are presented for: 

a. Case I in Figure 11 

b. Case II in Figure 22 

3. The optimal control is depicted for: 

a. Case IA in Figure 12 

b. Case IB in Figure 15 

c. Case IC in Figure 18 

d. Case IIA in Figure 23 

e. Case IIB in Figure 26 

f. Case IIC in Figure 29 

4. The optimal and desired altitude trajectories are 
presented for: 



a. 


Case 


IA in Figure 13 


b. 


Case 


IB in Figure 16 


c . 


Case 


IC in Figure 19 


d. 


Case 


IIA in Figure 24 


e. 


Case 


IIB in Figure 27 


f. 


Case 


IIC in Figure 30 
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5. The components of the optimal trajectory, v, 
Ot, 0, and 0, together with the related airplane parameters 
rate of descent and normal acceleration, are presented for 

a. Case IA in Figure 14 

b. Case IB in Figure 17 

c. Case IC in Figure 20 

d. Case IIA in Figure 25 

e. Case IIB in Figure 28 

f. Case IIC in Figure 31 

In all cases, the results were within established 
specification limits and were considered realistic for an 
airplane landing. The anticipatory nature of the optimal 
control was evident in the results, especially in Case IIA 
where the optimal control produced a finite thrust pertur- 
bation at the initial time, although the airplane was as- 
sumed to be in an ideal flight condition. This apparent 
lack of continuity, along with the discontinuity encounter 
ed in the optimal control of Cases IB, IC, IIB, and IIC at 
t = 0, was predictable, since it was assumed that the con- 
trols were instantaneously available. Development of a 
state model to include control dynamics, presented in 
Chapter III, will presumably eliminate this apparent dis- 
crepancy. 
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PROBLEM AREA 



c . 

The importance of proper selection of the constant 
values for the elements of the diagonal weighting matrices 
of the performance measure is evident; namely, no unique 
values exist and hence proper selection is based on pro- 
ducing an admissible control and trajectory when the op- 
timization technique is applied. In this investigation, 
where many specifications had to be simultaneously satis- 
fied, this selection became paramount, especially when it 
became apparent that several components of the state tra- 
jectory were sensitive to changes in one element of the 
weighting matrices. As a result, the trial-and-error pro- 
cess required to obtain values for these weighting matri- 
ces was tedious and time-consuming. Further study is 
warranted toward the development of a scheme or method to 
aid the control engineer in the selection of these weight- 
ing matrices. 
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FIGURE 10 

THE SOLUTION OF THE 
RICCATI EQUATION FOR CASE I 
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FIGURE 12 
THE OPTIMAL 
CONTROL FOR CASE IA 
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time-; 

FIGURE 13 

THE OPTIMAL AND DESIRED 
ALTITUDE TRAJECTORIES FOR CASE IA 




71 



FIGURE 14 

THE OPTIMAL a, 9 , AND 8 
HISTORIES AND h , F TRACES FOR CASE I A 
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FIGURE 15 
THE OPTIMAL 
CONTROL FOR CASE IB 
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THE OPTIMAL AND DESIRED 
ALTITUDE TRAJECTORIES FOR CASE IB 
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FIGURE 17 

THE OPTIMAL a , 6 , AND Q 
HISTORIES AND n , ft TRACES FOR CASE IB 
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THE OPTIMAL a, 9 , AND 9 
HISTORIES AND n , ft TRACES FOR CASE IC 







o 

to 



2 ; 



CD 



CO 



CO 



LO 



ro 

- cvt 

I 

I 

- o 



78 



FIGURE 21 

THE SOLUTION OF THE 
RICCAT! EQUATION FOR CASE II 
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FIGURE 23 
THE OPTIMAL 
CONTROL FOR CASE II A 
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THE OPTIMAL AND DESIRED 
ALTITUDE TRAJECTORIES FOR CASE II A 
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time-sec 

FIGURE 25 

THE OPTIMAL v , a , 9 , AND 9 
HISTORIES AND n , R TRACES FOR CASE HA 
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FIGURE 26 
THE OPTIMAL 
CONTROL FOR CASE IIB 
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FIGURE 28 

THE OPTIMAL v , a , 9 , AND 9 
HISTORIES AND n . fi TRACES FOR CASE I IB 
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FIGURE 29 
THE OPTIMAL 
CONTROL FOR CASE II C 




FIGURE 30 

THE OPTIMAL AND DESIRED 
ALTITUDE TRAJECTORIES FOR CASE IIC 
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FIGURE 31 

THE OPTIMAL v , a , 6 , AND 9 
HISTORIES AND n , F> TRACES FOR CASE IIC 



VIII. CONCLUSIONS 



For the results presented in Chapter VII, the follow- 
ing conclusions are stated: 

1. The feasibility of obtaining realistic results by 
formulating the all-weather landing problem as a linear 
tracking problem and applying optimal control techniques 
was demonstrated. 

2. Satisfactory landings, using a simple mathemati- 
cal model of the F-4J airplane, were accomplished in all 
cases investigated when the derived optimal control was 
applied. 

3. Obtaining values for the elements of the weighting 
matrices, H, Q, and R, to produce admissible and realistic 
results was difficult because the state trajectory was 
sensitive to changes in these matrices. The process was 
problematical, since a large number of design specifications 
had to be satisfied. 

Recommendations for logical extensions of this investi- 
gation follow: 

1. Formulate the problem to include the control dy- 
namics. 

2. Formulate the problem to include measurement noise 
and wind effects. 

3. Investigate application of sub-optimal control 
techniques to the problem. 
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4. Investigate the possibility of using the sensi- 
tivity of the state trajectory to changes in the weighting 
matrices as a basis for developing a method of obtaining 
representative values for the elements of these matrices. 
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APPENDIX A REPRESENTATIVE DATA FOR THE F-4J AIRPLANE 



IN THE LANDING CONFIGURATION 



Gross Weight 


32,500 lb 


Equilibrium Velocity (v ) 


223 ft/sec 


Stall Velocity (V g( _ all ) 


195 ft/sec 


Equilibrium Angle of Attack (a Q ) 


0.33 rad 


Stall Angle of Attack (£ gta n) 


0.49 rad 


Limit Sink Rate at Touchdown 


23 ft/sec 


Approximate Thrust Available 


MIL engine operation 


5,000 lb 


A/B engine operation 


12,000 lb 


Aerodynamic Ground Effect 


Negligible 


Control Time Constants 


Elevator (Tg ) 

P 


0.10 sec 


W 

Thrust (Tg ) 

e 


0.15 sec 


Dimensional Stability Derivatives 


X v =-0.593xl0 -1 sec 1 


Z & = °*° 


X = 0.0 


Zg =-0.13 2x10 2 ft sec -2 


q 


°e 


X a = 0.107xl0 2 ft sec 2 


Zg = 0 . 289xl0 -3 f t lb -1 


“n 

o 

• 

o 


-3 -1 

M v = 0.450x10 ft sec 


*6= °-0 


M =-0.430 sec 1 
q 


Xg =0.950xl0 -3 ft lb -1 sec -2 
1 


-2 

=-0.130x10 sec 


Z v =-0.259 sec -1 


=-0.260 sec 1 


Z =0.0 

q 


Mg e =-0.279 sec -2 


7 -7 


Z =-0*827x10 ft sec 

OL 


M* = 0.0 
°T 
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APPENDIX B SUMMARY OF NUMERICAL VALUES 



STATE MODEL CONSTANT 



a ll 


= -0.59300x10 1 


a 52 


= 


-0.22300x10 


a l2 


= 0.10700xl0 2 


a 53 


= 


0.22300x10 


a l3 


= -0.32172x10 


b 12 


— 


0.94667x10 


a 21 


= -0 . 11628xl0 -2 


b 21 


= 


-0.59193x10 


a 22 


= -0.37085 


b 22 


= 


0.12979x10 


a 41 


= 0.75232xl0~ 3 


b 41 


= 


-0.27746x10 


a 42 


= -0.12036x10 


b 42 


= 


-0.33745x10 


a 44 


= -0.69000 


C 51 


= 


-0.11676x10 


DIAGONAL WEIGHTING MATRIX 


CONSTANTS 






CASE I 










h ll = 1 *° 


q ll 


= 


1.0x10 1 




h 22 = 1,0 


q 22 


= 


1 . OxlO -1 




h 33 = 1.0 


q 3 3 


= 


1.0x10 1 




h 44 = 2 .OxlO -3 


q 44 


= 


5. OxlO -4 




R =5.0 










CASE II 










h-^ = 5.0x10 3 


q ll 


= 


1 . OxlO -5 




h 22 = 5. OxlO -1 


q 22 


= 


1.0x10 1 




h 33 = 5. OxlO -1 


q 33 


= 


1.0x10 1 




h 44 = 1.0 


q 44 


= 


5. OxlO -1 




h 55 = 5. OxlO -3 


q 55 


= 


5. OxlO -4 




r ll = 5 *° 


r 22 


= 


5. OxlO -10 



3 

3 

-3 

-1 

-5 

-6 

2 
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COMPUTER PROGRAM 



DERIVATION OF THE OPTIMAL CONTROL 
FOP AN ALL-WEATHER AIRPLANE LANDING SYSTEM 

NPGS MASTER THESIS, JUNE 1969 
MAIN PROGRAM 



REAL *8 ITITLEI 12) ,KKDCT(20) » KK ( 20 ) , XX DOT (5) » XX ( 5 ) • 
*ITITL1<36) ,nTK,TOO.TFF,DTX 
REAL *4 A< 25 ) ,B< 12) ,H(5) ,0(5) ,R (4) , RRTF<5) ,AT (25 ), 
*K(2C) ,S( 5) ,RRT(5) .TEMPI (25) .TEMP2 ( 25 ) , TEMP3 ( 25 ) , 

+QC( 25 ) , SDCT ( 5 ) , KDOT ( 25 > . X ( 5 ) , G ( 5 ) ,F(15).CC(5).P<15>. 
*SVALS( 5,1001 ) , XVALS <5,i601 ) ,USTAR(2, 1001 ) ,TBA$E( 1001) , 
♦BT(12),RI (A) , LAB (30) , HH C 5 ) ,HHH(5) .TEMP4(25) ,L1,L2,L3, 
♦GBASE ( 600 ) ,ALFDCT( 1001) ,ALTDOT( 1001 ), ACCEL ( 1001) , 
*KVALS( 15 ,1001 ) ,XXX< 1 00 1 ) , PSC AL ( 600 ) 

GRAPHICAL OUTPUT DATA. 

ITITLl(I) ARE GRAPH TI/TLES 
LAB ( I ) ARE CURVE LABELS 

R EAD (5,524) ( I T I TLl ( I ) ♦ I =1 , 36 ) 

WP I TE ( 6, 505 ) 

WRITE(6,525)(ITITL1(I ),I=1.36) 

READ( 5,531) ( LAB ( I ) , I =1 , 1 5 ) 

READ (5,531) ( L Ae ( I ) , I = 16 . 30 I 

WRITE! 6,507) ( L AB ( I ) ,1=1, 15) 

WRITE(6,5C7) (LAB( I ) ,1=16,30) 

INITIALIZING DATA 

NCASE IS CASE NUMBER ( 1 = 1 .2 = 1 1 ) 

NO IS NUMBER OF PLANT STATES 

NC IS NUMBER OF CONTROLS 

NSTK IS NUMBER OF INTERGRATION STEPS 

DT IS INTEGRATION STEP SIZE 

TO IS INITIAL TIME 

TF IS FINAL TIME 

L1,L2,L3 ARE CONSTANTS OF DESIREC TRAJECTORY ECUATION 



R EAD ( 5 , 5C C ) 
WRITE (6, 501) 

P E AD ( 5 , 5C2 ) 
DTX=DT 
DTK=-DTX 
WP I TE ( 6 , 5C3 ) 
EPS=0. IE-04 
NT=C 

NK= N0*( NO+1 )/2 
NKT=NK+NC 
NONO= NO*NC 
NCNC = NO*NC 
NSTK 1 = NSTK + l 
NSTK2 =NSTK+2 
NSTK 3=NS TK /2 
NSTK4=NSTK3+1 
NSTK5=NSTK+3 
R AD= 57.2958 
GR AV=3 2. 1725 
GAMMA= -3.0/RAD 
VC= 223.0 
ALTO=100.C 



NCASE .NG.NC, NSTK 
NCASE, NO, NC, NSTK 
DT,T0,TF,L1,L2,L3 



DT,T0,TF,L1 ,L2,L3 



TBASE IS TIME 
TIME BASE FOR 



BASE FCR COMPUTATION WHILE GBASE IS 
GRAPHICAL OUTPUT 



TBASE( 1 ) =0 . 0 
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OC 100 I =1 »NSTK 
ICO TEASE ( I+1)=TBASE( I ) + DT 
DC 101 J = 1 , NSTK4 
101 GPASE < J ) = -TBASE ( 2*J-1 ) 

PRGBLEM DATA 

A ( I ) ARE ELEMENTS CF A MATRIX 
B(I) ARE ELEMENTS CF B MATRIX 



,25) 

, 10 ) 
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READ ( 5 , 5C4 ) (A(I),I=1 
READ ( 5, 504 ) ( B 1 1 ) , I = 1 
DC 102 1=1, NO 

102 CC ( I )= 0.0 
CC(NO)=VC*GAMMA 
IF (NCASE.EO. 2 ) GO TO 
DO 103 1=1, NO 
8 ( I ) = B ( I ♦ 1 ) 

A ( I ) =A( I +6 ) 

103 A ( 1 + 4 ) = A ( I+li) 

DC 104 1=1, NO 
A ( I+8)=A( 1+16) 

104 A( I+12)=A( 1+21 ) 

105 WRITE ( 6, 508 ) 

DO 106 J=1 ,NQ 

106 WRITE(6,50S)(A(I) ,I=J, 

WRITE( 6 , 510 ) 

DO 107 J=1 , NO 

107 WRITE (6,509) (B( I) , I =J ,NCNC , NO ) 
WRITE (6, 511) 

WRITE ( 6,516 ) (CC(I), 

R E AD ( 5 , 504 ) (X(I),I = 1 
WRITE(6,523> 

WRITE (6, 516) (X(I),I=] 

RE AD ( 5 , 504 ) (H(I),I=1 
REAC ( 5 , 504 ) (Q(I),I=1 
READ ( 5 , 5C4 ) ( R ( I ) • I = 1 
RE AD ( 5 , 504 ) (RRTF(I), 
WPITE(6,512) 

WRITE! 6,505 ) <H(I),I=] 

WR ITE ( 6, 509 ) (0(1), 1=] 
WRITE(6,5C9) ( R ( I ) , I = 1 
WR ITE ( 6 , 522 ) 

WRITE (6,516 ) ( RRTF ( I ) , 1 = 1 ,N0 ) 



I , NCNQ , NO ) 



» 1=1, NO) 
. , NO ) 

! 1 , NO ) 

, NO ) 
l , NO ) 

. , NO ) 

J = 1 ,N0) 

a, no 
a, no 
a , no 



FORM A TRANSPOSE, B TRANSPOSE, AND R INVERSE 

CALL MTRA ( A , AT , NC , NO , 0 ) 

CALL MTPA(B,BT ,N0,NC,0) 

IF(NC.EO.l) GO TO 109 
CALL MSTR(P, TEMPI, NC, 2,1) 

CALL SINV (TEMPI, NC, EPS, IER) 

CALL MSTR(TEMP1,RI,NC,1,2) 

IF ( IER- 1 ) 110,108, 108 
109 WRITE ( 6, 514 ) 

GC TO 300 

109 R I ( 1 )= 1.0/Rdl 

110 WRI TE ( 6, 526 ) 

DO 111 J=1 , NO 

111 WRITE (6, 509) ( AT( I ) , 1= J , NONO , NO) 

WRITE ( 6, 527) 

DO 112 J=1 , NC 

112 WRITE (6, 509) ( BT ( I ) , 1= J ,NONC , NC ) 

WRI TE ( 6 , 528 ) 

WRITE! 6,509) ( R I ( I ) ,I=1,NC) 



FORM FINAL TIME VALUES FOR K AND S 



DO 113 1=1, NO 
113 RRTF ( I ) =-RRTF ( I ) 

CALL MSTR(H,K,NC,2,1) 

CALL MPRD(H,RRTF,S»N0,N0»2»0,1) 
WRITE(6,515) 
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WRITE (6,516) ( K( I ) ,I=1,NK) 

WRITE(6,517) 

WP I TE ( 6 , 5 16 ) ( S( I ) ,I=1,N0) 

DO 114 1=1 » NO 

114 RRT ( I ) =-PRTF ( I ) 

FORM P MATRIX WHICH EQUALS B*Rl*BT 

CALL MPRD(B,RI ,TEMP1,NC,NC,0,2,NC) 

CALL MPPD(TEMPl,BT,TEMP2,N0,NC,0,0,NO) 
CALL MSTR(TEMP2»P ,N0 » 0 » 1 I 
WRITE(6,513) 

WPITE(6,516) (P( I ) ,1=1, NK) 

CALL MSTR(0»QQ,N0,2»0) 

FORM KOOT AND SDCT AND INTEGRATE BACKWARDS 

DO 115 1=1, NK 

115 KK ( I ) =K ( I ) 

DO 116 1=1, NO 

116 KK ( NK-f I )=S( I ) 

Of 131 L= 1 , NSTK 



THE DESIRED 
TF, DT.NSTK 
FOLLOWING 9 



TRAJECTORY 
OR OESIREO 
CARDS MUST 



IS INTRODUCED HERE. I 
TRAJECORY ARE CHANGEC, 
BE CHANGEC ACCORDINGLY 



117 

118 
ns 
120 

121 

122 



123 



124 



125 

126 



127 

128 

129 

130 

131 

132 



IF(L.GE.401) GC TO 117 
DELT= TBASE(401)-TBASE(L) 

TEMP= EXP(L1*DELT)-1.0 

RR T ( NO ) = ( 1.0/Ll) *L2*TEMF - L2*CELT 

PRT(NO-l)=L3*DELT 

RRT (NO-2 )=(L3/2.0)*DELT**2 

R R T( NO-3 )=RRT(NC-2)-( L2/VO)*TEMP 

IF (NCASE.EQ. 1 ) GO TO 119 

RRT ( NO-4 ) = 0.0 

GC TO 119 

00 118 1=1, NO 

RP T ( I ) = 0.0 

DO 120 1=1, NK 

KVALSt I ,L)=K< I ) 

DO 121 1=1, NO 
SVALS ( I , L ) = S( I ) 

CALL MPRD( K »A, TEMPI, NC, NO, 1,0, NO 
CALL MPRD ( AT,K,TEMP2,NG,NO,0,1 ,N0) 

CALL MPRD (K,P,TEMP3,NC,NO,l,l.NO) 

CALL MPRD( TEMP3,K,TEMP4,N0,N0,0,l,NO) 

DC 123 1=1 » NONO 

KDOT ( I )= -TEMPI (I )-TEMP2( I > — 0 0 C I ) +TEMP4 ( I ) 
CALL MS TR (KDOT, TEMPI, NC,0,1) 

DO 124 1=1, NK 
KKDOT( I )= TEMPI ( I ) 

CALL MSUB (AT, TEMP3.TEMP4.N0, NO, 0,0) 

CALL MPRD(TEMP4,S,TEMPl,N0,N0,0,0,l) 

CALL MPRD(Q,RRT,TEMP2,Nu,Nu,2,0,l) 

DO 125 1=1, NO 

SOOT ( I )=-TEMPl( I ) +TEMP2 ( I ) 

DO 126 1=1, NO 

KKOOT ( NK+I ) = SDGT(I) 

SS=RKLDEQ(NKT,KK,KKOOT,TFF,DTK,NT) 

IF ( SS- 1. ) 127,128,128 
WRITE( 6,521 ) 

GO TO 300 
DO 129 1=1. NK 
K( I )= KK ( I ) 

00 130 1=1, NO 
S( I )= KK ( NK+ I ) 

IF ( SS- 1 . ) 127,122,131 

CONTINUE 

DO 132 1=1 , NK 

KVAL S( I.NSTK + l ) = K ( I ) 

DO 133 1=1, NO 



F TO, 
THE 
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133 S VAL S ( I ,NSTK+i )=S< I ) 

WRITE K AND S VALUES 

N I =1 

IF(NK-IO) 134,134,135 

134 NR=NK 

NT I = 1 
NR I =G 
GC TO 136 

135 NTI=NK/10 
NR=1C 

NR I =NK- NTI *1 0 

136 DO 138 M=1 ,NTI 
WR ITE ( 6 ,5 1 8 ) 

00 137 J=1,NSTK1 

137 WR I TE ( 6,509 ) ( K VALS ( I ,NSTK2-J)»I=NI,NR) 
NI=NI+10 

138 NR=NF+10 

1 FI NR I ) 141.141,139 

139 WRITE!6,518) 

NPI=NRI+NTI*10 
DO 140 J=1 ,NSTK 1 

140 WR ITE ( 6, 509 ) ! KVALS! I ,NSTK2-J ) , I = NI , NRI) 
IC=0 



10 READ (5.10C0,END=999) CARD 
1000 FORMAT I20A4I 



141 WRITE! 6,52C I 
DO 142 J= 1,NSTK1 

142 WRITE! 6,509 ) ( SVALS ( I »NSTK2-J)» 1 = 1, NO) 

OBTAIN OPTIMAL CONTROL AND OPTIMAL TRAJECTORY 



NT=0 

DO 143 1=1, NO 

143 XX! I )=X( I ) 

XXDOT ( NO-3 ) = 0.0 

XXDOT ( NO I = VO*GAMMA+VC* ( X ( NO-2 )-X l NO-3 ) ) 

DC 144 1=1, NO 

144 HH ( I )=0.0 

DO 157 L=1 , NSTK 
DO 145 1=1, NO 

145 XVALS! I ,L»=X! I J 

ALTDOT IS THE RATE OF DESCENT AND ALFDOT IS THE 
DERIVATIVE OF ALPHA 

ALTDOT1L)=XXDOT!NO» 

ALFDOT! L > = XXDOT ( NO-3 » 

147 DC 148 1=1, NK 

148 TEMP 1 ( I )= KVALS! I ,NSTK2-L» 

DC 149 1=1. NC 

149 TEMP2! I )= SVAL S I I .NSTK2-L ) 

CALL MPRD! PI ,0T,TEMP3.NC,NC.2,O,NC) 

CALL MPRD 1 TtMP3, TEMPI ,F , NC , NO.O , l , NC ) 

CALL MPRD 1 TEMP3,TEMP2,G,NC,NC, 0,0,1) 

HH(NO) IS EQUILIBRIUM ALTITUDE TRAJECTORY AND 
HHH IS F*HH 

HH! NO I=ALTO -( VC*3. C/RAO ) *TBASE ( L ) 

CALL MPRD! F , HH.HHH , NC , NO , 0,0,1 ) 

CALL MPRD! F , X , T EMP3 , NC , NO ,0 ,0 , U 
DO 150 1=1, NC 

15C TEMP4! I )= -TEMP3! I)-G< I ) ♦HHH! I ) 

DO 151 1=1, NC 

151 USTAR ! I , L ) = TEMP4 ! I I 

152 CALL MPRD! A, X, TEMPI, NC, NO, 0,0,1 ) 

CALL MPRD (B, TEMP4,tEMP2, NO, NC, 0,0,1) 

DO 153 I =1 ,N0 

153 XXDOT 1 I ) = TEMPI ( I )♦ TEMP2 ! I ) +CC ( I ) 
SS=RKLDEQ(NO,XX,XXDOT,TOO,OTX,NT) 
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IFC SS — 1* ) 154,155,155 
1 54 WP 1 TE ( 6 • 52 1 ) 

GC TO 300 

155 DO 156 1=1. NO 

156 X ( I ) = XX( I ) 

I F ( SS- 1 • ) 154,152,157 

157 CONTINUF 

DC 156 1=1, NO 

158 XVALS (I.NSTK+1) =X(I) 

AL TOOT ( N$TK + 1 ) *XXD0T(N0) 

A L FOOT (NSTK+l)=XXD0T(NC-3) 

WRITE OPTIMAL CONTROL AND OPTIMAL TRAJECTORY 

16C WPITE( 6,529) 

DO 161 J=1 ,NSTK 

161 WRITE! 6,530) TBASE ( J ) , ( USTAR ( I , J ) , I* 1 , NC ) 

ACCEL IS THE NORMAL ACCELERATION 

DO 162 J=1,NSTK1 

162 ACCEL! J)=VO*(XVALS(NO-l,J)-ALFDCT(J) J/GRAV *1.0 
WRITE! 6,519) 

00 163 J=1,NSTK1 

163 WRITE! 6,530) TBASE ! J ) , ( XVALS < I , J ) , I = 1 ,NC ) , ALTOOT ! J ) , 
♦ACCEL ( J ) 

GRAPHICAL OUTPUT SUBROUTINE CRAW , DEVELOPED 

AT NPGS » IS USED TO CONSTRUCT THE FOLLOWING PLOTS 
(SUBROUTINE AVAILABLE AT NPGS COMPUTER FACILITY) 

1. PLOT OF OPTIMAL CONTROL 

2. PLOT OF OPTIMAL AND DESIRED ALTITUDE 
TRAJECTORIES 

3. PLOT OF VELOCITY, ALPHA, THETA, THETA DOT, 
NORMAL ACCELERATION AND RATE OF DESCENT 

4. SOLUTION TO RICATTI EQUATION 

5. PLOT OF S VALUES 

DC 1 1=1,12 

1 ITITLE! I )=ITITL1( I ) 

DO 7 I = 1 , NC 
IF(NC.FO.l) GO TO 4 
GO TO! 2 , 3 ) ,! 

2 MCD=1 
GF S = 1 • 0 
GO TO 5 

3 MCD=3 

GP S=C. 0GC1 
GO TO 5 

4 MOD=C 
GR S= 1. 0 

5 DC 6 L =1 , NSTK3 

6 XXX! I )=USTAR( I ,2*L-1)*GRS 

7 CALL DRAW(NSTK3, XXX, GBASE, MOD. 0, LAB(I ), ITITLE, 

*0 .05, l.C, 11, 4, 2, 2, 9, 11,1, LAST) 

DC 8 1=7.12 

8 ITITLE! I » = I T I TL 1 ( 1+6) 

DO 15 1=1,2 

GC TO (9, 13), I 

9 MOD= 1 

IF DESIRED ALTITUDE TRAJECTORY IS CHANGED, NEXT 
SIX CARDS MUST PE CHANGEC ACCORDINGLY 

DC 10 L = 1 ,600 

10 XXX ( L ) = ALTO+VO*GAMMA*TBASE(L) 

DO 11 L=6C1 »NSTK1 
DELT=TBASE(L)-TEASE(601) 

TEMP=EXP!L1*DELT>-1.0 

11 XXX! L)=ALTO+VO*GAMMA*TBASE(L ) +L2*TEMP/L1-L2*DELT 
DO 12 L=1 , NSTK4 

12 XXX! L ) = XXX! 2*L-1 ) 

GO TO 15 

13 MGD=3 
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00 14 L = 1 , NSTK4 

14 XXX(L) = XVA LS ( NO » 2*L-1 ) 

15 CALL DRAW( NSTK4 ,XXX ,G8ASE ,MOD ,0 ,LAB(I*2) , ITITLE, 

*20. 0,1. 0,1 1,1, 2, 2, 9, 11,1, LAST) 

DO 16 1=7,12 

16 ITITLE (11= ITITLHI + 12) 

DO 27 1=1,6 

GO TO ( 17,18,18,18,23,25) , I 

17 MOD= 1 
JJ=0 

IF ( NCASE. EQ. 1 ) GO TO 19 
GRS=0.0 1 
GO TO 21 

18 M0D=2 
GRS=1.0 
GO TO 21 

19 DO 20 L=1,N$TK4 

20 XXX( L ) =0. 0 
JJ = 1 

GO TO 27 

21 DO 22 L= 1 » NSTK4 

22 XXX( L ) = XVALS( I- J J ,2*L-1 ) *GRS 
GO TO 27 

23 DC 24 L=1 , NSTK4 

24 XXX ( L ) = ACCEL(2*L-1)*0.1 
GC TO 27 

25 MCD=3 

DC 26 L=1 ,NSTK4 

26 XXX( L ) = ALTDOT(2*L-1)*0.01 

27 CALL DRAW( NSTK4, XXX, GBASE , MOD, 0 ,LAB( 1+4) , ITITLE, 

*0.05, 1.0,11 ,4,2,2,9,11,1 , LAST) 

DO 28 1=7,12 

28 ITITLE< I )=ITITL1( 1+18) 

DC 29 1 = 1 » NSTK4 

29 GBASE( I ) = -TBASE(NSTK5-2*I ) 

DC 37 1=1, NK 

DC 30 J=1 , NSTK4 

30 PSCAL ( J )= KVAL S ( I ,2*J-1 ) 

CALL SCALE(NSTK4. PSCAL, 60. 0,GRS) 

GO TO (31,32,32,32,32,32,32,32,32,33,32,32,32,32,34), I 

31 MOD= 1 

GO TO 35 

32 M0D=2 

GC TO 35 

33 I F ( NCASE. E 0. 1 ) GO TO 34 
GO TO 35 

34 M0D=3 

35 DO 36 L=1 , KSTK4 

36 XXX(L)=KVALS( I ,2*L-1)*GRS 

37 CALL DRAW(NSTK4, XXX, GBASE, MOD, 0 ,LAB( 1+10) , ITITLE, 
*60.0,1.0,11,4,2,2,9,11 ,1 ,LAST) 

DO 38 1=7,12 

38 I T I TLE ( I )= I TITL1 ( 1+24) 

DO 46 1=1, NO 

DC 39 J=1 ,NSTK4 

39 PSCAL( J )= SVAL S ( I ,2*J-1 ) 

CALL SC ALE(NSTK4, PSCAL, 3.0, GRS) 

GC TO (40,41,41,42,43),! 

40 MCD=1 

GO TO 44 

41 M0D=2 

GC TO 4<> 

42 I F ( NCAS F • EO. 1 ) GO TO 43 
GO TO 44 

43 M0D=3 

44 DO *5 L = 1 , NSTK4 

45 XXX( L )= SVALS( I ,2*L-1 )*GRS 

46 CALL DRAWINSTK4, XXX. GBASE, MOD, 0,LAB(I+25), ITITLE, 

*3. 0,1. 0,1 1,4, 2, 2, 9, 11,1, LAST) 

300 CONTINUE 

500 FORMAT (4110) 

501 FORMAT ( // / , T8 , • NCASE= • , 14, / , T8 , • NO* • , 14 , / , T8 , 
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* • NC = ' »I 4 , / , T8 . ’NSTK=’ ,14) 

502 FORMAT ( 6E 10. 5 I 

50 3 FORMAT (T8, • 0T= • , E 13 . 5 , / , T8 . 'T0= • , El 3. 5 , / , T8 , *TF=* , 
♦E13. 5,/,T8,’Ll=’,E13.5,/,T8,'L2=',El3.5,/,T8, 

*’TF=* ,F13.5) 

504 FORMAT ( 5E 10. 5 ) 

505 FORMAT ( ’1’ ,//) 

506 FORMAT ( 8E 1 3* 5 ) 

507 FORMAT(2X, 151A4.2X) ) 

508 FORMAT(///»T9, 'A MATRIX’,/) 

509 FORMAT! 10E13.5) 

510 FORMAT! ///,T8, ’6 MATRIX’,/) 

511 FOP MAT! ///,T8,’C MATRIX’,/) 

512 FORMAT!///, T8. 'H,Q, AND R MATR IC ES ( D I AG STORAGE)’,/) 

513 FORMAT !///,T8,'P MATRIX (SYMETRIC STORAGE)’*/) 

514 FORMAT!///, T8, ’ERROR IN SINV’,/1 

515 FORMAT!///, T8,’K(TF) ! SYMETR I C STORAGE)',/) 

516 FORMAT(E 13. 5) 

517 FORMAT!///, T8. 'S(TF) VECTOR’,/) 

518 FORMAT ( ///, T8 , *R ICATTI COEF',/) 

519 FORMAT ( ///,T8 ,’ STATE SOLUTIONS',/) 

520 FORMAT <///,T8,’ S COEF', /) 

521 FORMAT ! ///, T8» ’ERROR IN RKLDEQ ' ,/ ) 

522 FORMAT! ///,T8, ’RR(TF) VECTOR’,/) 

523 FORMAT!///. T8, ’X(O) VECTOR’,/) 

524 FORMAT ( 6A8 ) 

525 FORMAT ( 2 X ,{ A8 ) 

526 FORMAT! ///,T8, ’A TRANSPOSE ’,/ ) 

527 FORMAT! ///,T8, *6 TRANSPOSE’,/) 

529 FORMAT!///, T8, »R INVERSE ( DI AG STORAGE)',/) 

529 FORMAT!///, T8, ’OPTIMAL CONTROL’,/) 

530 FORMAT(F6. 2.8E14.5) 

531 FORMAT ( 1 5A^ ) 

END 



FUNCTION RKLDEQ 



PURPOSE 

THIS ROUTINE SOLVES A SYSTEM OF N FIRST-ORDER 
ORDINARY DIFFERENTIAL EQUATIONS BY THE RUNGE-KUTTA- 
GILL FOURTH-ORDER METHOD 

DESCRIPTION OF PARAMETERS 

N - NUMBER OF FIRST-CROER EQUATIONS IN SYSTEM TO 
BE SOLVED. (O.LE.N.LE. 25) 

Y - NAME CF LINEAR ARRAY OF LENGTH AT LEAST N, 

IN WHICH SCLUTION VALUES WILL 8E STOREO. THE 
CALLING PROGRAM SHOULD SUPPLY INITIAL VALUES 
BEFORE FIRST ENTRY. 

F - NAME CF LINEAR ARRAY OF LENGTH AT LEAST N, 

IN WHICH THE DERIVATIVES, COMPUTED IN USER'S 
PROGRAM APE STOREC. 

X - THE INDEPENDENT VARIABLE, WHICH IS ACVANCED 
WITHIN! RKLOEQ 

H - THF INCREMENT FOR X, WHICH MAY EE CHANGED AT 
THE END OF ANY INTERVAL. ! WHEN S=2.0) 
t! - AN INTEGER WHICH COUNTS THE NUMBER CF TIMES 
ENTRY TO RKLDEQ HAS BEEN MAOE DURING CURRENT 
INTERVAL. IT MUST BE I NIT ALLY SET TC ZERO BY 
USER BtFORE FIRST CALL OF RKLOEQ. SUBSEQUENTLY 
IT SHOULD NOT BE CHANGED BY USER. 

S - A SWITCH TC BE TESTED BY USER UPON RETURN FROM 
RKLDEQ. 

IF S = 1.0, THE CALLING PROGRAM SHOULD NOW 
COMPUTE VALUES OF F, USING CURRENT VALUES OF 
X AND Y, AND RETURN TO RKLDEQ. 

IF S - 2.0, THE END OF PRESENT INTERVAL HAS 
BEEN REACHEO. USER SHOULO STORE AND/OR OUT- 
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PUT CURRENT X AND/CR Y AND TEST FOR END OF 
COMPUTATION. 



FUNCTION RKLDEC ( N , Y , F , X , H, NT ) 

REAL*8 Y»F»X,H,G»H1»H2»H3»H6 
DIMENSION Y ( 1 1 , F(l), 0(25) 

NT = NT +1 

GC TO (1,2, 3, 4), NT 

1 HI = H 

H2 = HI * • 5D0 
H3 = HI * 2. DO 
H6 = H1/6.D0 
DO 11 J =1 , N 
11 O(J) = 0. DC 
A = . 5D0 
X = X + H2 
GO TO 5 

2 A = .2928932188134525 
GO TO 5 

3 A = 1.7071C67811865475 
X = X + H2 

GO TO 5 

4 DO 41 I = 1,N 

41 Y ( I ) = Y ( I ) + H6 * F ( I I -Q<II/3.00 
NT = 0 
RKLDEO =2. 

GO TO F 

5 DC 51 L = 1 , N 

Y(L> = Y ( L ) ♦ A * ( H * F ( L ) -0(D) 

51 0 ( L ) = H3 * A *F(L) ♦ < 1 . DO-3 . DO * A ) *0(U 
RKLDEO =1. 

6 RETURN 
END 



SUBROUTINE SCALE 



SUBROUTINE SCALE ( NPTS ,PTS , SCLE ,GRS ) 
DIMENSION PT S ( 600 ) 
RANGE=3.0*SCLE+SCLE/4.0 
MA X=ABS ( PTS ( 1 ) ) 

DC 11 I =2, NPTS 
TNAX=ABS( PT S ( I ) ) 

IF(MAX-TMAX)10,ll,ll 
1C MA X= AB S ( P T S ( I ) ) 

11 CONTINUE 

I F ( RANGE-MAX I 12,22,13 

12 KFLAG=1 

IF(MAX.LE.(RANGE*10.0) ) GO TO 14 
IF(MAX.LE. (RANGE*100.0) ) GO TO 15 
IFIMAX.LE. (RANGE*1000.0) ) GO TO 16 
IF(MAX.LE. (RANGE*1C000.0) ) GO TO 17 
IF(MAX.LE. (RANGE*100000.0)) GO TO 18 
GC TO 18 

13 KFLAG=2 

IF(MAX.GE. (RANGE*0. 1) ) GO TO 22 
IF(MAX.GE. (RANGE*0.01 ) I GC TO 14 
IF(MAX.GE. (RANGE* 0.001) ) GO TO 15 
IF(MAX.GF. (RANGE*0.0001 ) ) GO TO 16 
IF(MAX.GE. (RANGE*0. 00001) ) GO TO 17 
IF(MAX.GE. (RANGE*0. 000001 )) GO TO 18 
GO TO 16 

14 TGRS=10. C 
GO TO 1° 

15 TGRS=100. 0 
GO TO 19 

16 TGRS=100C. 0 
GC TO 19 

17 TGRS=10000.0 
GO TO 19 
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10 TGRS= 100000. 0 
19 GO TO < 20,21 ) , KFLAG 
2C GPS= 1.0/TGRS 
GC TO 23 

21 GRS= TGRS 
GO TC 23 

22 GR S = 1 • 0 

23 CONTINUF 
RETURN 
END 
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